function [ W, H ] = CNMF( duomenys, R, niter, varargin)

    M = size(duomenys, 1);
    N = size(duomenys, 2);
    
    Y = duomenys;
    W0 = rand(M, R);
    H0 = rand(R, N);
    
    [alpha beta W0 H0] = Popt(varargin, 'alpha', 1, 'beta', 0, 'W0', W0, ...
                                    'H0', H0);
    W = W0;
    H = H0;
    
    ee = 1e-6;
    
    for k=1:niter
        W = W .* (Y*H')./((W*H)*H' + alpha*W + ee);
        H = H .* (W'*Y )./(W'*W*H + beta*H + ee);
    end

end

